Dynamical Properties of Two Coupled Hubbard Chains 

at Half-filling 
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Using grand canonical Quantum Monte Carlo (QMC) simulations combined with Maximum En- 
tropy analytic continuation, as well as analytical methods, we examine the one- and two-particle 
dynamical properties of the Hubbard model on two coupled chains at half-filling. The one-particle 
spectral weight function, A(k, u), undergoes a qualitative change with interchain hopping t± asso- 
ciated with a transition from a four-band insulator to a two-band insulator. A simple analytical 
model based on the propagation of exact rung singlet states gives a good description of the features 
at large t±. For smaller t±, A(k,u;) is similar to that of the one-dimensional model, with a coherent 
band of width the effective antiferromagnetic exchange J reasonably well-described by renormal- 
ized spin-wave theory. The coherent band rides on a broad background of width several times the 
parallel hopping integral t, an incoherent structure similar to that found in calculations on both 
the one- and two-dimensional models. We also present QMC results for the two-particle spin and 
charge excitation spectra, and relate their behavior to the rung singlet picture for large t± and to 
the results of spin-wave theory for small t±. 



I. INTRODUCTION 

The compounds SrCu 2 O 3 and (VO) 2 P 2 7 i consist of 
arrays of weakly interacting two-chain metal-oxide lad- 
ders. In SrCu203, Cu0 2 ladders are coupled by a weak, 
frustrated, ferromagnetic couplinga, and in (VO) 2 P 2 07, 
VO4 ladders are well separated in the structure of the 
material. These materials have a half-filled conduc- 
tion band, corresponding to one conduction electron per 
metal-oxide unit and are insulators with short-range an- 
tiferromagnetic correlations and a spin gap. Very re- 
cently Z. Hiroi and M. Takanda succeeded in doping the 
two-leg ladder compound Lai_ 2; Sr 2 ;Cu0 2 .5. They ob- 
serve a fall in resistivity with doping, possibly to a metal- 
lic state, but no sign of superconductivity. The spin- 
gap found in half-filled LaCu0 2 .5 persists, but decreases 
with doping. In this publication we will address the 
half-filled situation only. The undoped materials have 
been described by a Heisenberg model on two coupled 
chainstffl "13, which explains a number of experimentally 
observed magnetic properties of the materials, including 
the size of the spin gap, and the spectrum of the spin 
triplet excitations. The Hubbard model on two coupled 
chains, which at half-filling and in the limit of strong on- 
site Coulomb repulsion can be mapped onto the Heisen- 
berg model, gives an itinerant electron picture of these 
systemsB which includes charge as well as spin fluctua- 
tions. While the ground state and the low-energy proper- 
ties of the two-chain Heisenberg and half-filled Hubbard 
model are becoming well understood, less is known about 
the finite frequency one- and two-particle response. The 
dynamic spin response of the Heisenberg model was stud- 
ied in Ref. 5 and the single particle and spin response was 



examined for a t—J model with two holes in Ref. 8 using 
analytic continuation of Lanczos exact diagonalization 
calculations on small clusters. The dynamic properties 
are important because they can be measured with a va- 
riety of experimental techniques. For example, the one- 
particle spectral weight can be probed by photoemission 
and inverse photoemission, the dynamic charge correla- 
tion function by optical response measurements, and the 
dynamic spin structure factor by inelastic neutron scat- 
tering. 

Here we examine the dynamic one- and two-particle 
correlation functions of the two-chain Hubbard model 
at half-filling using Quantum Monte Carlo simulations 
analytically continued tojjeal frequencies using a Max- 
imum Entropy techniqueE3. We study the dynamic re- 
sponse for large, intermediate, and small values of the 
perpendicular hopping t±. This is interesting for two 
reasons: first, when the electron-electron interaction is 
turned off (U — 0) in the two-chain system, there is 
a metal-insulator transition from a two-band metal to 
a band insulator due to the separation of the bonding 
and antibonding bands with increasing t±. For the in- 
teracting (U 7^ 0) case, it is generally believed that the 
two-chain system is always insulating (see, for example, 
Ref. 11 for a discussion within a weak-coupling RG pic- 
ture). The QMC simulations do find an insulator for 
all t±, but also reveal a crossover from four-band insu- 
lating behavior for small t± to two-band behavior for 
t±/t > 2. A similar crossover occurs for weak coupling 
(U t) within antiferromagnetic H artre e-Fock (AFHF) 
theory, as we shall discuss in section [II B| . One interesting 
and important result of this paper is that this four-band 
to two-band crossover scenario survives in the interme- 
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diate to large U regime in the system, despite the fact 
that the many-body physics for t ± > t is not even qual- 
itatively reproduced in the AFHF approximation. No 
such crossover is present as a function of J±/J in the 
Heisenberg model, because the mapping to the Heisen- 
berg model breaks down when J± ~ &b\/U is of the 
order of t±. In other words, the effect of the electronic 
band structure is not present in the Heisenberg model. 
As we shall see, this crossover can also be described in 
terms of the crossover from a localized rung picture, in 
which localized excitations form a Bloch state along the 
chains for large t±, to a picture with longer range anti- 
ferromagnetic correlations with many features that can 
be qualitatively described by AFHF spin-wave theory for 
small t±. 

Secondly, in the small to intermediate ij_ regime, this 
system displays features of the single-particle spectral 
weight A(k,ui\ seen in rpcpnt numerical treatments of 
both the lDBB and 2lMt5 Hubbard model. We find 
two features also generic to the ID and 2D systems: a 
dispersive band whose width is of the order of the effec- 
tive exchange interaction J ~ 4t 2 /U and an incoherent 
background of width several times the parallel hopping 
integral t. For the one and two chain systems, but not 
the 2D system, the dispersive band is well described by 
renormalized spin-wave-theory. 

We consider the single band Hubbard model on two 
coupled chains of length L: 

H =- t Y,( c iM c i+iM +h - c -) 

z, Act 

-t± ^(ct lo .c ii2(r + h.c.) + n iM n iM - (1) 

Here c 1 - > and c- , create and annihilate electrons on 

l, ACT I. ACT 

rung i and chain A with spin a, respectively, the hopping 
integral parallel to the chains is t, the hopping between 
the chains t±, and U is an on-site Coulomb interaction. 
We use periodic boundary conditions parallel to and open 
boundary conditions perpendicular to the chains. 

The non-interacting, U — 0, Hamiltonian can be 
diagonalized by writing the hopping term in terms of 
bonding and antibonding states on a rung and Fourier- 
transforming parallel to the chains. The energy is then 
given by 

£k = — (2icos k + t± cos fcjj (2) 

with k = (fc, k±), where fcj_ = and fcj_ = n corresponds 
to the energy of the bonding and antibonding band, re- 
spectively, and k is the momentum along the chains. 
Both bands will be occupied when tx < ij_ c , whereas 
when ti_ > <j_c only the bonding band will be occupied. 
Here t± c /t = 1 - cos7r(n) and (n) = (J2a 4,x.a c ^A,cr) ■ 
For half-filling t± c /t = 2, and the system is a two-band 
metal for tx/t < 2, and a band insulator with a com- 
pletely occupied bonding band for tx/t > 2. 



II. SINGLE-PARTICLE SPECTRAL WEIGHT 

A. Quantum Monte Carlo Results 

In order to calculate the dynamical properties of the 
above model we employ the grand-canonical Quantum 
Monte Carlo (QMC) algorithm to determine the expec- 
tation values of the correlation functions in imaginary 
time. At half-filling, the simulations are not limited by 
the fermion sign problem and accurate results can thus 
be achieved at low temperatures and large system sizes. 
We analytically continue the imaginary-time data to real 
frequencies using a Maximum Entropy methodEl For ex- 
ample, the spectral weight function given by 

l,v 

■ 8{w- {E v - Ei)), (3) 

where Z is the partition function, \l > and \V > are 
the exact many-body eigenstates with the corresponding 
energies Ei and E v and c k , CT = J2i.x c i ,\cre lk '' R *- x / V%L, 
can be calculated by inverting the spectral theorem 

. r°° e -™ 

G(k,r) = (c k|T (r)ci tT (0))= / AQt, U )dw 

(4) 

with the Maximum Entropy algorithm using the raw 
QMC data G(k, r). To achieve high resolution, it is 
important to use a likelihood function which takes the 
error-covariance matrix of the QMC data-and its statis- 
tical inaccuracy consistently into accountllj. The results 
presented here are based on QMC data with good statis- 
tics, i.e. averages over 10 5 updates of all the Hubbard- 
Stratonovich variables result in G(k, r)'s with absolute 
errors less than or of the order of 5 x 10~ 4 . Correlations 
of the data in imaginary time were taken into account 
by making use of tha covariance matrix in the Maximum 
Entropy; procedurdlS. As suggested in previous work by 
WhiteEZl, various moments of the spectral weight were 
also incorporated in extracting A(\s.,u>). 

We calculate the spectral function A(\s., u>) at half- 
filling for a 2 x 16-lattice with U/t = 8 for different val- 
ues of t± at an inverse temperature of — 10/i. Our 
results are given in Fig. [II for t±/t = 2.0, in Fig. |2| for 
tx/t = 1.0, and in Fig. | for t ± /t = 0.5. Note that 
tx/t = 1.0 corresponds to the isotropic case, J± w J, 
relevant to the SrCu203, (VO)2P207, and LaCu02.5 
compounds. Parts (a) and (b) in Fig. are three- 
dimensional plots of A(k, u>) versus u> for fc-values in the 
ID Brillouin zone, whereas part (c) in all three figures 
summarize these results in the usual "band structure" 
uj versus k plot. The shaded areas represent regions 
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with significant spectral weight, with the darker shad- 
ing representing more weight. Due to particle-hole sym- 
metry, A(k = (k,ir),uj) = A(k = (tt — k,0),— lo) where 
k = (k, k±). In other words, one reflects k about k = tt/2 
and lo about to get A(k, lo) for k± = tt from A(k, lo) 
for k±_ = 0. This symmetry can be seen by comparing 
parts (a) and (b) of Fig. Since this symmetry is not 
enforced by the Maximum Entropy procedure, it is only 
present approximately, and the extent to which the re- 
flected spectra at k± = do not match those at k± = tt 
gives an indication of the accuracy of the analytic con- 
tinuation procedure. In the density plots in parts (c), 
we show only the k± — components. For t±/t = 2.0, 
there is a heavily weighted coherent band of width ~ it 
in the photoemission (lo < 0) part of the spectrum. In 
the inverse photoemission spectrum (u> > 0) there is very 
little total spectral weight. There are therefore two bands 
with significant spectral weight, one in the photoemission 
spectrum for k± = and one in the inverse photoemis- 
sion spectrum for k± — tt. Each band is about 2t away 
from the Fermi surface, leading to a gap of approximately 
At. Therefore, for large t±, the system is a two-band in- 
sulator. 

In contrast, for t±/t = 1.0 and t±/t = 0.5, A(k, lo) 
has substantial spectral weight in four bands, one in the 
photoemission spectrum (lo < 0) for k±_ = and k± = tt 
and one in the inverse photoemission spectrum (lo > 0) 
for both values of k±. In this regime, the system is thus 
a four-band insulator. For t±/t = 0.5 (Fig. ||) the spec- 
tral weight is present for both values of k± concentrated 
between k = and k = tt/2 in the photoemission part 
and between k = tt/2 and k = tt in the inverse pho- 
toemission part of the spectrum. For the isotropic case 
(t±/t = 1.0, Fig. ||), there is some shift of spectral weight 
to the photoemission part for k± = and to the inverse 
photoemission part for k± = tt. The maxima of the pho- 
toemission band in the k± = part occurs at k* = tt for 
tj_/t = 2.0, at k* w 0.7?r for t±/t = 1.0, and at k* = tt/2 
for t±/t — 0.5. Therefore we would expect to see a max- 
ima at k* w Q.Itt in a photoemission experiment on the 
isotropic ladder compounds, an experiment which has, to 
our knowledge, not yet been done. 

It is also important to note that the QMC A(k, lo) 
results for t±/t = 0.5 as well as for the isotropic case, 
t±/t = 1.0, display two general features which have re- 
cently been observed in both the lDt2 and 2DEil cases 
using the impxpved Maximum Entropy techniques de- 
scribed aboveEj. One feature is that A(k, lo) contains 
a rather dispersionless "incoherent background" extend- 
ing over several t (~ 6t for U/t = 8 in 2D) in both the 
electronically occupied (lo < 0) and unoccupied (lo > 0) 
parts of the spectrum. The crucial structure, not re- 
solved in previous ID and 2D QMC simulations, is a 
dispersive structure at low energies with a small width 
of the order of the exchange coupling J = At 2 /U . It 
is this latter coherent "band" which defines the gap A 
and which, for larger U (U/t ~ 10), is well-separated 
from the higher-energy background. As in the ID and 



2D cases, the splitting in the low-energy band and the 
higher-energy background is especially pronounced near 
k = and k — tt due to a relative weight shift from 
negative to positive energies as k moves through fc*. 

B. Spin— Wave Theory 

At half-filling, the two chain Hubbard system at T = 
is a spin liquid with a spin gap and a finite spin-spin 
correlation length. The system is therefore less ordered 
than either the ID system, which is at the critical point 
and has power law decay of the spin-spin correlation 
function and no spin gap, or the 2D system, which has 
long-range spin order and gapless spin-wave-like excita- 
tions. However, at large U/t the scale of the spin gap 
for the two chain system is set by J± ~ At 2 j_/U, and 
the spin-spin correlation length £ grows longer as the 
spin gap gets smaller. The range of the spin ordering 
can therefore be tuned by varying t±/t. This is illus- 
trated by Fig. U in which we plot the spin-spin corre- 
lation function, S(r) = (— l) r (Mp A M r 2 A ) on a semilog 
scale for tj_/t = 2.0, tj_/t = 1.0 and ii/i = 0.5. Here 
Af^ A = (?v.A,f — fV,A,i) is twice the z component of the 
on-site spin. These calculations were done usi&g the Den- 
sity Matrix Renormalization Group methods (DMRG) 
on a system with open boundaries at zero temperature. 
For t±/t = 0.5, the data was averaged over a number 
of different pairs of points for a given r in order to re- 
duce oscillations due to the open boundaries. For dis- 
tances greater than a few lattice spacings, the correla- 
tion functions are very well fit by a pure exponential 
form Acxp(— r/£), from which the spin-spin correlation 
length £ can be determined. For t±/t = 2.0, £ = 0.83, 
for t±/t = 1.0, £ = 4.3, and for t±/t = 0.5, £ = 9.5 
lattice spacings. Therefore, the spins are almost com- 
pletely uncorrelated along the chains for t±/t = 2.0, but 
for t±/t = 1.0 and for t±/t = 0.5 are correlated on length 
scales of the order of the size of the 2 x 16 lattice stud- 
ied here, although the long-range behavior is that of a 
gapped spin liquid in all three cases. 

It is therefore reasonable to calculate the single par- 
ticle spectral weight A(k, lo) in two simple ways: first, 
when the spin gap is small and £ « L, and we consider 
wavelengths smaller than £ and frequencies larger than 
the spin gap, it is useful to explore the consequences of 
a simple SDW approximation based on an ordered state, 
i.e. AFHF theory. Secondly, one can solve for the ex- 
act eigenstates on a rung and then consider states com- 
posed of a product of exact rung singlets only and a de- 
localized one-hole state in a background of rung singlets. 
(We call this the Local Rung Approximation, LRA.) The 
LRA starts with a state with uncorrelated rungs and then 
treats the interactions between rungs perturbatively, and 
is thus a good starting point when £ is small. As we shall 
see, the LRA gives a single particle spectral weight distri- 
bution that agrees well with QMC in the large t± regime, 
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while AFHF qualitatively describes most features of the 
spectral weight distribution for small t±. 

We first discuss the AFHF calculation in more detail. 
The Hartree-Fock one-particle Hamiltonian can be writ- 
ten as 

^f = E'^K>m-CC). (5) 

k.o- 

where means that the sum runs only over the mag- 

netic zone (— 7r/2 < k < tt/2, k± = 0, 7r). The operators 
a an< ^ a are gi ven by the following transformation!!!!]: 

b La = u kCw,a T VwCW+Qm, (6) 
K,<r = VkC k:!T ± UkCk+Q : <r, (7) 

where Q = (tt, tt), the upper (lower) sign corresponds to 
o =T (c =i), and 



«w = y 5 (i+ek/^), (8) 



v k =<\/-Q.-e k /E k ), (9) 

£ k = VA 2 + £ k 2 . (10) 

Here £k is the free particle energy ({/ = 0) defined in Eq. 
(||). The spin-density-wave gap, A, is self-consistently 
determined by the equation 

k 

The band structure of the AFHF Hamiltonian is given 
by ±£k = ±\/A 2 + £k 2 in the magnetic Brillouin zone. 
The SDW gap A, calculated by solving Eq. (0) self- 
consistently, is shown as a function of t±/t for U/t = 2, 
4, and 8 in Fig. The transition from a two-band metal 
to a band insulator in the U = system occurs at tj./t = 
t± c = 2.0. As U is turned on for t < t± c in the AFHF, a 
gap develops in both the bonding and antibonding bands, 
leading to a four-band insulator. If t± is then increased, 
A will go to zero at approximately the noninteracting 
(U = 0) tj_ c /t = 2.0, as seen for U/t = 2 in Fig. |, 
leading to a transition from a four-band to a two-band 
insulator. For large U and small t±, A w C7/2, as seen 
for f7/f = 8 in Fig. ||, so the Coulomb splitting of the 
bands will dominate over {7 = band structure effects. 
In this case, the gap will go to zero, i.e. the system 
will undergo a transition to a band insulator only when 
the band structure splitting 2t±, is of the order of the 
Coulomb splitting, so that t± c s» A w [7/2. 

The spectral weight within the AFHF approximation 

is 

A(\£,w) = ul5(w-E ] J+vl5(L> + Ek). (12) 

We show the spectral weight A(k, oj) from AFHF for the 
same parameters as in Fig. |l| and Fig. || in Fig. ^(a) for 



t±/t = 2.0 and Fig. |(b) for t±/t = 0.5. The band split- 
ting 2A sa U in both cases, as can also be seen in Fig. ||. 
Since both the bonding and antibonding bands are split 
by this gap, there are two peaks in A(k, uS) as a func- 
tion of u> for each k, leading to a four-band insulator for 
both values of t±/t. Since the band splitting is set by U 
rather than 2t± , the spectral weight is almost evenly dis- 
tributed between two coherent bands for both t±/t= 2.0 
and t±/t — 0.5. Therefore, for t±/t = 2.0, as can be seen 
by comparing the QMC results in Fig. [l] with Fig. ||(a), 
the AFHF spectral weight distribution for t±/t = 2.0 is 
completely wrong, with too much weight in the upper, 
lo > 0, band for k± = and the uj < band for k± = tt. 
While the average positions of the AFHF bands, shown as 
thick solid lines in Fig. 0(c) , are approximately the same 
as those of the QMC bands, the band widths are too 
small by approximately a factor of two. For t±/t = 2.0 
and U/t = 8, then, the AFHF calculation results in a 
four-band insulator in which there is long-range anti- 
ferromagnetic order, while the QMC results show that 
the system is a two-band insulator, with only very short 
range antiferromagnetic correlations. While the AFHF 
picture does produce a two-band insulator for t±/t > 
4.5, the physical picture of the transition to this phase is 
quite different (see section C). 

For t ± /t = 0.5, as seen in Fig. | and Fig. |(b), AFHF 
gives two dispersive bands in the k±_ = branch of width 
of order J ~ 4t 2 /U, similar to the coherent bands seen 
in the QMC data. However, in the AFHF, the single 
particle gap is somewhat too large, the weight distribu- 
tion extends too far towards k = tt and too far towards 
k = for the photoemission and inverse photoemission 
parts of the spectrum, respectively, and the broad inco- 
herent background is not present. We iave also carried 
out mean-field slave boson calculationsEJ in order to try 
to improve on the AFHF picture. These calculations give 
a gap about 15% smaller than AFHF, and a bandwidth 
about 3% smaller, but do not qualitatively improve on 
the AFHF calculations. The same holds for t±/t = 1.0, 
where the AFHF describes approximately the coherent 
part of the spectrum and also produces a too large single 
particle gap and a slightly different spectral weight dis- 
tribution. Therefore, for t±/t = 0.5 and for t±/t = 1.0, 
AFHF would give a reasonable description of the disper- 
sion and general spectral weight distribution of the coher- 
ent spin-wave bands, if the gap were phenomenologically 
adjusted to a smaller value, but fails to even qualitatively 
describe the spectral weight distribution for t±/t = 2.0. 

C. Local Rung Approximation 

For very large t±, a better starting point is the limit 
of weak interaction between the rungs. In this limit, we 
split the Hamiltonian of Eq. (Q) into H = Hq + Hj with 

i.er iX 
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where Hq denotes the non-interacting rung limit. The 
ground state of Hq is a product of individual rung eigen- 
states. Accordingly, we first diagonalize the Hubbard 
Hamiltonian on the two sites making up the rung in or- 
der to find the rung eigenstates. A schematic diagram of 
the exact eigenstates for a single rung is shown in Fig. 
0. The ground state for the half-filled rung i (two par- 
ticles per rung) is a spin singlet state with k± = and 

energy E a /L = E™ ng = (u - y/U 2 + 16^) /2 denoted 

by \Si). For large values of U, the energy of this two-site 
state is given by the exchange coupling — Jj_ ~ —At\/U. 
Note that \Si) contains terms which have double oc- 
cupied sites as well as the usual singlet construction 

(IU)-U,T»/V2. 

We can then form an approximation to the ground 
state for the half-filled system by taking a state which is 
just a product of the rung states, 

\^) = \S x )\S 2 )...\Sl). (14) 

For t±/t = 2.0, the binding energy of a singlet formed 
on a rung is about four times lower than the energy of 
a singlet between neighboring sites along a chain, and 
therefore we expect |^o) to be a good approximation to 
the exact ground state. 

In order to obtain the u < spectral weight for the 
half-filled system, we need to calculate the matrix ele- 
ment of Ck,o- given in Eq. (0). In other words, we need 
matrix elements between a half-filled state and a state 
with one particle removed. We form the approximate one 
hole state by replacing one bond singlet state at rung I, 
\Sg), with the lowest energy one-particle state, the one 
with bonding (k± = 0) symmetry, \Bi^) (We remove a 
spin down electron for definiteness.) and define 

\£) = \S 1 )\S 2 )...\B £] )...\S L ). (15) 

We delocalize the single-particle r s±ate with a plane wave 
ansatz by constructing the state&EII 

L 

hMfc)) = £- 1/2 E e ^>- ( 16 ) 

Trivially, this state is an exact eigenstate of Hq. We 
will take this state as the ground state for the system 
with one hole with momentum k = (k, 0). The spectral 
weight .A(k, ui) at zero temperature (/3 — > oo) is then ap- 
proximated by using only the two states |^q) and \ipi(k)} 
in Eq. (0). For k± = 0, the energy dispersion u>(k) of 
A(k, ui) for lu < is given by 

u(k) = (^o|#|Vo) - (Mk)\H\Mk)) - H 

= -- ^U 2 + 16t 2 ± +t±-tA cos k, (17) 



w ith A = ( 1 + E 2 /2t ± ) 2 /{l + E 2 /At\), E 2 = (U + 
\J\J 2 + 16tj~)/2, and the corresponding spectral weight 
by KV'i(fc)l c fe,j.lV ; o)| 2 - The dispersion given by Eq. ( |T7j ) 
for U/t = 8 and t ± /t = 2.0 is plotted in Fig. 0(c) and 
the spectral weight for k±_ — is shown in FigTgl. One 
can see that the position of the peak of the u> < LRA 
band, denoted LRA1 in Fig. 0(c), lays almost exactly on 
the QMC data. However, the position and the disper- 
sion of the lu > band, which were obtained in the same 
way as for to < 0, are not exactly matched by the LRA1 
calculation. While this band is not very important in 
the sense that it contains very little spectral weight, we 
can understand how to improve the LRAl calculation by 
considering the states of a four site cluster. 

In order to generate the inverse photoemission spec- 
trum, we need to calculate matrix elements between the 
half filled state and a state with one additional particle. 
From Fig. 0, we see that a rung state with three parti- 
cles and fc_L = has approximately twice the excitation 
energy of the k± = one-particle state which is relevant 
for the photoemission part of the spectrum. (Recall that 
due to the particle -hole symmetry at half- filling, this 
state will map to one in the inverse photoemission part 
of spectrum for k = it; one can see this symmetry in Fig. 
0.) Since the relevant three-particle single rung state is 
high in energy, configurations involving intrachain effects 
might also be important. We include the effect of such 
configurations by replacing two rung singlet states by the 
lowest energy state of five particles on four sites in the 
bonding channel. We can then form a delocalizcd plane 
wave state from this state as in Eq. (|l6|). The results of 
this calculation are labeled as LRA2 in Fig. 0(c). The 
location and width of the k± — inverse photoemission 
band are more closely fit than by the LRAl calculation. 
In Fig. H the spectral weight distribution of the LRAl 
calculation for the photoemission part of the spectrum 
(u) < 0) and the distribution of the LRA2 calculation for 
the inverse photoemission part (lu > 0) is plotted. The 
result is in good accordance with the QMC data in Fig. 
0(a). 

III. SPIN AND CHARGE DYNAMIC 
CORRELATION FUNCTIONS 

In order to determine the nature of the low lying exci- 
tations, we also consider the spin and charge susceptibil- 
ities Xs,c(q, uj) which, in a Lehmann representation, are 
defined as 

Xa,c(q,«) = ^E e ^ Si ( 1 - e-^)|(/|O s , c (q)|0| 2 
M' 

• S(u - {E v - Ei)), (18) 

with 0,(q) = E P (cJ,+q,T c p,T - cj, +q x c p x ) and O c (q) = 
J2 P a c p+ q .cr c p,cr- We calculate the two-particle dynamic 
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response from the two-particle imaginary time Green's 
function, as in Eq. (Q) and use the Maximum Entropy 
method described above for the analytical continuation 
of Xs,c(q><£) to real frequencies. 

In Fig. we show the dynamical spin-spin correlation 
function Xs,c(<l> w) for t±/t = 2.0 and t±/t = 1.0. For 
t±/t= 2.0, the q±_ = component [Fig.||(a)] has abroad, 
lightly weighted structure centered approximately at lj = 
3t, with the spectral weight vanishing for small q. For 
q± = tt [Fig. |](b)], there is a coherent, dispersive band 
with a width set by J « At 2 /U = 0.5t, a minimum at 
q = tt and a maximum at 5 = 0. The minimum spin gap 
is approximately 0.8t, which agrees well with the value we 
obtain using Projector Quantum Monte Carlo, indicated 
by the dashed line on the plot. The spectral weight is 
most heavily concentrated around q = it. 

The dynamic spin response at half-filling measures the 
response of the system to spin triplet excitations. Since 
here t± is large, we consider the effect of triplet excita- 
tions on product rung states of the type considered in 
Eq. dbf), in the LRA calculation. A triplet excitation 
will change the total spin of the state from S = Q to 
5 = 1. Referring to Fig. fj], the only triplet excited state 
on a single rung has momentum it (corresponding to odd 
parity under chain interchange) leading to a change in 
momentum qj_ = tt for the triplet excitation. The size 
of this triplet excitation A££ ung is marked on Fig. ||(b) 
by the solid line. This local triplet excited state can be 
moved to a near neighbor rung by a process which is sec- 
ond order in Hj. The local triplet excited states can then 
be delocalized into a Bloch wave as in Eq. (|l|) and Eq. 
(H). For large U, one then obtains a dispersion rela- 
tion of the form AE r s ung + Jcosq, where J ~ At 2 /U , a 
form similar to that obtained in Ref. 5 for the two chain 
Heisenberg model. The coherent band in Fig. ^(b) does 
have a minimum at q = it, and has approximately this 
form. 

Within the LRA picture, it is also possible to excite 
a q± = (even under chain interchange) triplet excita- 
tion via a "two-magnon" process, as discussed in Ref. 
5. This corresponds to making two local triplet excita- 
tions on rungs in the non-interacting rung picture, and 
would lead to an excitation energy whose scale, in the 
Heisenberg limit, is set by lj ~ 2J± rather than lj ~ Jj_. 
The overall position of the lightly weighted band seen for 
q±_ = in Fig. |^(a) is consistent with this two-magnon 
band, although the band is too lightly weighted and our 
resolution too low to extract a dispersion relation. 

For the physical relevant isotropic case (t±/t = 1.0), 
shown in Fig. |(c) and (d), the spin response looks quite 
different. (Because the t±/t = 0.5 results are qualita- 
tively similar to the t±/t =1.0 results, we show only the 
isotropic (tj_/t = 1.0) case here.) In Fig. ^|(c) and (d), 
there are dispersive bands in both the qj_ = and qj_ = tt 
branches, with the position of the peak going to a finite 
minimum at q = (0, tt) and q = (vr,0), and appears to 
vanish as q — * (tt, tt). However, as shown in Fig. [|, the 
correlation length is 4.3 lattice spacings in this regime, 



and there should be a spin gap. Using DMRG calcula- 
tions on lattices of 2 x 8 to 2 x 32 sitesQ, we estimate the 
spin gap to be of order 0.12£. Due to finite size effects 
and finite resolution in the analytic continuation, this 
gap is too small to resolve in Fig. ||(d). Near q = (0, 0) 
[Fig. §(c)], the dispersion of the peak is hard to discern 
because there is very little spectral weight. This lack of 
spectral weight at q = (0, 0) is present in all the dynamic 
charge and spin correlations and is due to a selection 
rule that comes about because s , c (q = 0) in Eq. ( |lg| ) 
commutes with the Hamiltonian, leading to a vanishing 
of the matrix element of O s . c (q) as q — > 0. Therefore, 
the system has relatively long-range spin order in this 
regime and, to within the resolution of our calculations, 
shows the characteristics of an ordered state with gapless 
excitations. Thus it is interesting to compare the QMC 
spectra to the results of spin-wave theory calculations. 

In order to extract the low-lying spin excitations us- 
ing spin-wave theory it is necessary to consider the RPA 
transverse spin susceptibility XrpaC 4 , k', lj). It is ob- 
tained by applying the random-phase approximation to 
the response function Xo~(k, k', lj), which is calculated 
directly in real time using the AFHF ground state |f2), 
calculated in section [IB: 



Xo + -(q,q',0 = ^(n|TS+(t)sz q/ (o)|n), (19) 

with S ± — S x ± iSy. Due to the broken spin rotational 
symmetry of Xrpa( c 1' l'> u ) contains a gapless mode, 
as predicted by the Goldstcpr: theorem. Following the 
procedure of Schrieffer et al.E£l, one obtains: 



Xrpa(<W> w ) = 

]Tx^(q,p,^)[i - tfx3"-(p,q»]-\ 



(20) 



where [1 — t/xo~ (p, q',w)] 1 is a matrix inverse of a 2 x 2 
matrix in q-space and 



Xa 



"(q,q',w) = <Kq-q')xo (q, 



-*(q-q' + Q)Xo~(q.w). 



(21) 



with Xo (q> w ) ano - Xq (q? w ) given by the usual "bub- 
ble" diagrams and printed in detail in Ref. 16 [up to a 
misprint of a factor of 2 in XQ _ (q: w )]- The dispersion 

of Xa,pA(q> q' u ) IS indicated as a solid line in Fig. ^|(c) 
and y(d). The RPA spin-wave dispersion goes to zero at 
q = (0, 0) and q = (tt, tt), and goes to a finite minimum 
at q = (0,7r) and q = (tt,0), consistent with the QMC 
data. The spin-wave velocity, the width of the bands, 
and the gaps at q = (0,7r) and q = (tt,0) are also in 
reasonable agreement with the QMC data. 

The QMC results for the dynamic charge correlation 
function are shown in Fig. 10 for t±/t = 2.0 and t±/t = 
1.0. For t±/t = 2.0, almost all of the spectral weight 
is in the q±_ = tt component, so we do not show the 
q± = component. There exist two important features: 
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one is a heavily weighted, relatively flat band at w ~ Wt 
with heaviest weight near q — 0. This band becomes 
somewhat incoherent as q increases. The second is a flat, 
dispersive, less heavily weighted band with a minimum 
of order 5t near q = n. The spectral weight in this band 
extends from about q = 7r/4 to q = tt, and the size of the 
charge gap is set by this lightly weighted lower band. 

In the large t± limit, one can understand the structure 
of the ch arge response from the LRA picture described 
in section II C. In the single rung picture, a charge exci- 
tation will occur through a transition to an excited state 
conserving the number of particles on the rung and the 
total spin; in other words, the important transition will 
be in the middle column of Fig. |?], from the low-lying 
5 = state to the higher 5 = states. There are two 
possible charge excited states on the rung, one with mo- 
mentum k± = (even parity) and one with momentum 
k± = 7r (odd parity). In a single rung picture, an optical 
transition from the k±_ = ground state to a k±_ = 
excited state is forbidden because the q± — density 
operator O c (0) commutes with the Hamiltonian. This 
selection rule forbidding a q± = optical transition re- 
mains present when the rung charge excited states are 
constructed as in Eq. (|l^) and delocalized in a state like 
that in Eq. (pf). This is why there is almost no spectral 
weight in the q± = portion of the charge response for 
t±/t = 2.0. Of course, the system is not exactly in a LRA 
state, so there will be some higher order processes that 
will introduce a very small amount of spectral weight into 
the q± = branch. The energy of the q± = ir single rung 
transition, indicated by a solid line on the Fig. fio|(a), 
gives an excitation energy that agrees well with energy 
of the heavily weighted region at q = 0. 

In order to qualitatively understand the origin of the 
dispersion of the charge response, one can consider the 
possible particle-hole excitations within the one-particle 
band structure given by A(k, oj) in Fig. |l|. There will be 
significant amplitude in the two-particle charge response 
when there is significant amplitude for a transition at a 
particular q = {q,q±) and ui for a particle-hole excita- 
tion built up from the one-particle spectral weight. For 
example, to understand the heavily weighted amplitude 
at q = (0,7r), one has to integrate over all transitions 
from the photoemission band in Fig. |l](a) to the inverse 
photoemission band in Fig. |l|(b) which transfer this mo- 
mentum. Since the single-particle bands are sharp and 
parallel, one should obtain a single well-defined peak for 
q = (0,7r), which we see in Fig. [to](a) by considering 
excitations between the u> < 0, k±_ = band, and the 
lu > 0, k± = 7r band. In addition, the transition at 
q = (0, 7r) has odd parity between the chains, and is thus 
allowed by the selection rules for the density operator as 
q — > 0. As the parallel component q is increased, one can 
see that there will be a continuum of excitation energies 
which gets wider as q increases. The minimum excitation 
energy (position of the lower band) from this construc- 
tion is shown in Fig. |l0|(a) by the line with solid dots. 
However, the excitation energy obtained is consistently 



smaller than the energy of the lowest heavily weighted 
band from QMC. Using the single particle bands to con- 
struct the two-particle excitations is equivalent to calcu- 
lating the charge response using the lowest order "bub- 
ble" diagram, but with exact single-particle propagators, 
neglecting all particle-hole interactions. The particle- 
hole interactions on the rung, which are included in the 
rung eigenenergies, thus raise-the charge excitation en- 
ergy by a substantial amountEa. 

For t±/t = 1.0, the charge response looks quite differ- 
ent. As shown in Fig. |l^(b) and (c), there is substan- 
tial spectral weight for both q_\_ = and qj_ = tt. For 
q — » (0,0), the density operator has even parity, causing 
the matrix element and thus the spectral weight to van- 
ish, whereas at q — > (0,7r), the density operator has odd 
parity so that optical transitions are allowed and there is 
spectral weight. In both channels, the size of the charge 
gap is approximately At, and there is a broad structure 
of width ~ It. For q± = 0, most of the spectral weight 
occurs as a dispersive peak whose energy increases with 
increasing q, whereas for q± = n two peaks seem to con- 
tribute to the spectral weight distribution. The peaks 
are not well-defined enough over a range of q to extract a 
dispersion, but the upper peak is heavily weighted near 
q = 7r, at ui ~ 9t. 

We can qualitatively understand the broad incoherent 
structure of the charge response by considering particle- 
hole excitation in the single particle A(k, ui) in Fig. ||. 
There are four dispersive bands and a broad background, 
so there should be weight in both the q± = and q± = 
7r branches of the charge response, and broad structure 
above a minimum excitation energy, which we see in Fig. 
[To|(c) and (d). From the single particle bands in Fig. ||, 
one can estimate the minimum particle-hole excitation 
energy to be ~ At. In Fig. [l0|(c) and (d), the spectral 
weight near this minimum excitation energy is suppressed 
due to the particle-hole vertex. We have also carried out 
a calculation of charge response X h °pa( c 1' i'j w ) within the 
antiferromagnetic RPA approximation described above, 
using the SDW dispersion in Eq. ([h]), and also find a 
relatively broad structure above the charge gap for both 
q± = and q± — tt. We plot the minimum excitation 
energy of Xrpa^ Q'; w ) m Fig- 0(b) and (c) as lines 
with solid dots. This line is located in the middle of the 
broad band in both plots. 

The spin and charge response functions for t±/t = 0.5 
show the same general features as in the isotropic (t± jt = 
1.0) case and are therefore not shown here. The spin 
susceptibility Xs(q, is also qualitatively identical to 
the RPA result Xkpa^ *!> w ) with a smaller spin velocity 
than in the t±/t — 1.0 case. The charge susceptibility 
Xc(q, oj) shows a clear dispersive band centered around 
the low-lying RPA excitations in the q± = channel, 
whose energy increases with increasing q. In the q± = tt 
channel, there is a broad structure of width ~ 8t with 
again two peaks in the spectral weight distribution. 
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IV. CONCLUSION 



In summary, the single and two-particle dynamical 
properties of the two chain Hubbard model at half-filling 
can be understood by starting from two limits: the limit 
of non-interacting rungs treated exactly, which gives a 
good starting point for the large t± case for which the 
spin-spin correlation length along the chains is less than 
a lattice spacing, and the limit of an antiferromagneti- 
cally ordered state, which gives a good starting point for 
the small t ± case when the spin-spin correlation length is 
large. The dynamical properties in the two regimes look 
quite different. For the large t± regime, the remnants of 
the level transitions of the two site system representing a 
single rung, suitably broadened into bands, can explain 
the major features of the single particle, spin and charge 
responses. In the small t± case, calculations based on an 
antifcrromagnetically ordered starting point, such as an- 
tiferromagnetic Hartree-Fock theory and spin-wave the- 
ory give a good qualitative picture of the coherent spin- 
wave part of the single particle spectral weight and of 
the two-particle spin dynamic correlation function. In 
addition, there is a broad incoherent band in the single 
particle spectral weight similar to that found in recent 
numerical work on the ID and 2D systems. 

We also have shown the single particle spectral weight 
and the spin and charge response functions for the phys- 
ical relevant, isotropic case (t±/t = 1.0). The results are 
qualitatively similar to these in the small t±_ region and 
can therefore be understood from an antiferromagnetic 
Hartree-Fock or spin-wave theory picture. The real lad- 
der compounds have approximately the same coupling 
strength parallel to and perpendicular to the chains, and 
it will therefore be interesting to compare these isotropic 
t±/t = 1.0 results with future experiments. 
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FIG. 1. QMC result for the spectral density function 
A(k, uj) for all possible values of k = (k, k±) on a 2 x 16 sys- 
tem with U/t = 8 and t±/t = 2.0 at an inverse temperature 
of f3 — 10 /t. Parts (a) and (b) are three-dimensional plots on 
the LJ-k plane for k± — and k± = n, respectively. Part (c) 
is a density plot for the k± — branch with darker shading 
corresponding to higher spectral weight, and the points with 
error bars showing the position of the peaks. The thick solid 
line indicates the SDW result within AFHF at T = 0, the 
dashed line is the calculated local rung approximation (LRA) 
dispersion for local excitations on one rung (LRA1) and the 
long-dashed line shows the spectrum obtained by local exci- 
tations on a four-site cluster (LRA2). 

FIG. 2. The single particle spectral weight A(k,u)) as in 
Fig. 1, with the same parameters except with an isotropic 
hopping, t±/t = 1.0. In (c), the solid line is the SDW result 
calculated within AFHF. 



FIG. 3. The single particle spectral weight A(k, ui) as in 
Fig. 1, with the same parameters except with t±/t — 0.5. In 
(c), the solid line is the SDW result calculated within AFHF. 

FIG. 4. The spin-spin correlation function S(r) for 
t± = 2t, t± = t and t± = 0.5t calculated on a 2 x 32 lat- 
tice with U/t — 8 and open boundary conditions using the 
DMRG. The solid, dotted and dashed lines are asymptotic fits 
to the data using the simple exponential form Aexp(— r/£). 

FIG. 5. The SDW gap calculated within the AFHF ap- 
proximation at T = 0, as a function of t±/t on a 2 x 100 grid 
for U/t = 2,4,8. 

FIG. 6. Spectral weight calculated within the AFHF ap- 
proximation for (a) t±/t — 2.0 and (b) t±/t = 0.5 with 
U/t = 8, shown for k± = and k± — tv. 

FIG. 7. Exact eigenstates for a two-site system repre- 
senting a disconnected rung, with U = 8 and t± = 2.0. The 
horizontal axis is the number of particles and the vertical axis 
shows the total energy ((H — fiN)), where the chemical po- 
tential fi — U/2 at half-filling. Each level is labeled with the 
total spin S and momentum k±. 

FIG. 8. Spectral weight using LRA1 calculation for the 
ui < and the LRA2 calculation for the ui > parts of the 
spectrum for k± = 0. Here t±/t = 2.0 and U/t = 8. 



FIG. 9. Dynamic spin response function Xs(q, ^) for (a) 
t±/t = 2.0, q± = 0, (b) t±/t = 2.0, q± = it, (c) t±/t = 1.0, 
q± = 0, and (d) t±/t = 1.0, q± = it. Again, U/t = 8 on a 
2 x 16 lattice. The weight of sizable structures is represented 
by the strength of shading in the shaded regions and the posi- 
tion of the maxima by the points with error bars. In (b), the 
solid line indicates the spin gap AE™ ng of the two-site system 
and the dashed line the spin gap AEs for the 2 x 16 ladder 
at T = 0, calculated with a projector QMC algorithm. In (c) 
and (d), the solid lines are Xrpa( < 1' 1> u ) calculated using the 
SDW approximation given by Eq. (19). 

FIG. 10. The dynamic charge response function Xc(q, w) 
for the same parameters as in Fig. 7. In (a), t±/t = 2.0 
and q± = n, the solid line represents the charge excitation 
energy on a single rung Ai5™ ng , and the line with solid dots 
represents the minimum excitation energy constructed from 
the single particle spectral weight in Fig. 1. In (b) and (c), 
t±/t — 1.0, and q± = and 7r, respectively. The lines with 
solid dots represent the minimum excitation energy of the 
RPA charge susceptibility Xr°pa( ( 1) w ) calculated with the 
method of Ref. 16. 
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Dynamical Properties of Two Coupled Hubbard Chains 
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Using grand canonical Quantum Monte Carlo (QMC) simulations combined with Maximum Entro- 
py analytic continuation, as well as analytical methods, we examine the one- and two-particle 
dynamical properties of the Hubbard model on two coupled chains at half-filling. The one-particle 
spectral weight function, j4(k, to), undergoes a qualitative change with interchain hopping t± asso- 
ciated with a transition from a four-band insulator to a two-band insulator. A simple analytical 
model based on the propagation of exact rung singlet states gives a good description of the features 
at large t±. For smaller t±, j4(k, to) is similar to that of the one-dimensional model, with a coherent 
band of width the effective antiferromagnetic exchange J reasonably well-described by renormali- 
zed spin-wave theory. The coherent band rides on a broad background of width several times the 
parallel hopping integral t, an incoherent structure similar to that found in calculations on both 
the one- and two-dimensional models. We also present QMC results for the two-particle spin and 
charge excitation spectra, and relate their behavior to the rung singlet picture for large t± and to 
the results of spin-wave theory for small t±. 



I. INTRODUCTION 

The compounds SrC^Os 1 and (VO)2P2C>7 2 consist of 
arrays of weakly interacting two-chain metal-oxide lad- 
ders. In SrCu2C>3, CuC>2 ladders are coupled by a weak, 
frustrated, ferromagnetic coupling 3 , and in (VO)2P2C>7, 
VO4 ladders are well separated in the structure of the 
material. These materials have a half-filled conducti- 
on band, corresponding to one conduction electron per 
metal-oxide unit and are insulators with short-range an- 
tiferromagnetic correlations and a spin gap. Very recent- 
ly Z. Hiroi and M. Takano 4 succeeded in doping the two- 
leg ladder compound Lai-^Sr^CuO^ 5. They observe a 
fall in resistivity with doping, possibly to a metallic state, 
but no sign of superconductivity. The spin-gap found in 
half-filled LaCuC>2 5 persists, but decreases with doping. 
In this publication we will address the half-filled situati- 
on only. The undoped materials have been described by 
a Heisenberg model on two coupled chains 3,5-7 , which 
explains a number of experimentally observed magnetic 
properties of the materials, including the size of the spin 
gap, and the spectrum of the spin triplet excitations. The 
Hubbard model on two coupled chains, which at half- 
filling and in the limit of strong on-site Coulomb repul- 
sion can be mapped onto the Heisenberg model, gives an 
itinerant electron picture of these systems 8 which inclu- 
des charge as well as spin fluctuations. While the ground 
state and the low-energy properties of the two-chain Hei- 
senberg and half-filled Hubbard model are becoming well 
understood, less is known about the finite frequency one- 
and two-particle response. The dynamic spin response 
of the Heisenberg model was studied in Ref. 5 and the 
single particle and spin response was examined for a t—J 



model with two holes in Ref. 8 using analytic continuati- 
on of Lanczos exact diagonalization calculations on small 
clusters. The dynamic properties are important because 
they can be measured with a variety of experimental tech- 
niques. For example, the one-particle spectral weight 
can be probed by photoemission and inverse photoemis- 
sion, the dynamic charge correlation function by optical 
response measurements, and the dynamic spin structure 
factor by inelastic neutron scattering. 

Here we examine the dynamic one- and two-particle 
correlation functions of the two-chain Hubbard model 
at half-filling using Quantum Monte Carlo simulations 
analytically continued to real frequencies using a Maxi- 
mum Entropy technique 10 . We study the dynamic re- 
sponse for large, intermediate, and small values of the 
perpendicular hopping tj_. This is interesting for two 
reasons: first, when the electron-electron interaction is 
turned off" (U = 0) in the two-chain system, there is a 
metal-insulator transition from a two-band metal to a 
band insulator due to the separation of the bonding and 
antibonding bands with increasing tj_. For the interac- 
ting (U 7^ 0) case, it is generally believed that the two- 
chain system is always insulating (see, for example, Ref. 
11 for a discussion within a weak-coupling RG picture). 
The QMC simulations do find an insulator for all 
but also reveal a crossover from four-band insulating be- 
havior for small t±_ to two-band behavior for t±/t ^ 2. 
A similar crossover occurs for weak coupling (U ^ t) 
within antiferromagnetic Hartree-Fock (AFHF) theory, 
as we shall discuss in section II B. One interesting and 
important result of this paper is that this four-band to 
two-band crossover scenario survives in the intermediate 
to large U regime in the system, despite the fact that the 
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many-body physics for t±_ > t is not even qualitatively 
reproduced in the AFHF approximation. No such cros- 
sover is present as a function of Jj_/J in the Heisenberg 
model, because the mapping to the Heisenberg model 
breaks down when Jj_ ~ At\/U is of the order of t]_. In 
other words, the effect of the electronic band structure is 
not present in the Heisenberg model. As we shall see, this 
crossover can also be described in terms of the crossover 
from a localized rung picture, in which localized excitati- 
ons form a Bloch state along the chains for large to a 
picture with longer range antiferromagnetic correlations 
with many features that can be qualitatively described 
by AFHF spin-wave theory for small 

Secondly, in the small to intermediate t±_ regime, this 
system displays features of the single-particle spectral 
weight j4(k,w) seen in recent numerical treatments of 
both the ID 12 ' 13 and 2D 14 ' 15 Hubbard model. We find 
two features also generic to the ID and 2D systems: a 
dispersive band whose width is of the order of the effec- 
tive exchange interaction J ~ At 2 /U and an incoherent 
background of width several times the parallel hopping 
integral t. For the one and two chain systems, but not 
the 2D system, the dispersive band is well described by 
renormalized spin-wave-theory. 

We consider the single band Hubbard model on two 
coupled chains of length L: 



H = -t'S^(c*-, a-,, , +h.c.) 

-t±. X!( c li<7 c *\2<7 + hx -) + u J2 ni ^ ni ' x ^ w 



Here c] „ and c • ^ „ create and annihilate electrons on 
rung i and chain A with spin a, respectively, the hopping 
integral parallel to the chains is t, the hopping between 
the chains and U is an on-site Coulomb interaction. 
We use periodic boundary conditions parallel to and open 
boundary conditions perpendicular to the chains. 

The non-interacting, [7 = 0, Hamiltonian can be dia- 
gonalized by writing the hopping term in terms of bon- 
ding and antibonding states on a rung and Fourier- 
transforming parallel to the chains. The energy is then 
given by 



-(2t cos k + t j_ cos 



(2) 



with k = (k, where = and = 7r corresponds 
to the energy of the bonding and antibonding band, re- 
spectively, and k is the momentum along the chains. 
Both bands will be occupied when t±_ < tj_ c , whereas 
when t±_ > tj_ c , only the bonding band will be occupied. 
Heretj_ c /t = l-cos7r(n) and (n) = (J2 a c ],\,<j c iA,<?) ■ For 
half-filling tj_ c /t = 2, and the system is a two-band me- 
tal for ti_/t < 2, and a band insulator with a completely 
occupied bonding band for tj_/t > 2. 
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1. QMC result for the spectral density function 
for all possible values of k = (k, k±) on a 2 x 16 sy- 
stem with U/t = 8 and t±/t = 2.0 at an inverse temperature 
of fi = 10/t. Parts (a) and (b) are three-dimensional plots on 
the Lo-k plane for k± = and k± = ir, respectively. Part (c) 
is a density plot for the k± = branch with darker shading 
corresponding to higher spectral weight, and the points with 
error bars showing the position of the peaks. The thick solid 
line indicates the SDW result within AFHF at T = 0, the 
dashed line is the calculated local rung approximation (LRA) 
dispersion for local excitations on one rung (LRA1) and the 
long-dashed line shows the spectrum obtained by local exci- 
tations on a four-site cluster (LRA2). 



II. SINGLE PARTICLE SPECTRAL WEIGHT 

A. Quantum Monte Carlo Results 

In order to calculate the dynamical properties of the 
above model we employ the grand-canonical Quantum 
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FIG. 2. The single particle spectral weight A(k, u;J as m 
Fig. 1, with the same parameters except with an isotropic 
hopping, t±/t = 1.0. In (c), the solid line is the SDW result 
calculated within AFHF. 

Monte Carlo (QMC) algorithm to determine the ex- 
pectation values of the correlation functions in imaginary 
time. At half-filling, the simulations are not limited by 
the fermion sign problem and accurate results can thus 
be achieved at low temperatures and large system sizes. 
We analytically continue the imaginary-time data to re- 
al frequencies using a Maximum Entropy method 10 . For 
example, the spectral weight function given by 



Z ^ 



-I3E, 



(l + e-^)|(/|c k)T |/'}r 



(3) 



where Z is the partition function, |/ > and \V > are the 
exact many-body eigenstates with the corresponding 




FIG. 3. The single particle spectral weight A(k, u;J as m 
Fig. 1, with the same parameters except with t±/t = 0.5. In 
(c), the solid line is the SDW result calculated within AFHF. 



energies E\ and E v and c k)(7 = J2i a c i,\<J e 
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can be calculated by inverting the spectral theorem 
G(k, r) = (c k T (r) C t T (0)} = / A(k, w)<h> 

(4) 

with the Maximum Entropy algorithm using the raw 
QMC data G(k,r). To achieve high resolution, it is 
important to use a likelihood function which takes the 
error-covariance matrix of the QMC data and its stati- 
stical inaccuracy consistently into account 16 . The results 
presented here are based on QMC data with good stati- 
stics, i.e. averages over 10 5 updates of all the Hubbard- 
Stratonovich variables result in G(k,r)'s with absolute 
errors less than or of the order of 5 x 10 -4 . Correlations 
of the data in imaginary time were taken into account 
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by making use of the covariance matrix in the Maximum 
Entropy procedure 16 . As suggested in previous work by 
White 17 , various moments of the spectral weight were 
also incorporated in extracting A(h,ui). 

We calculate the spectral function A(k,w) at half- 
filling for a 2 x 16-lattice with U ft = 8 for different va- 
lues of t± at an inverse temperature of /? = 10/t. Our 
results are given in Fig. 1 for t±/t = 2.0, in Fig. 2 for 
t L ft = 1.0, and in Fig. 3 for t±/t = 0.5. Note that 
ti_/t = 1.0 corresponds to the isotropic case, Jj_ J, 
relevant to the SrCu2C>3, (VO)2P2C>7, and LaCuC>2.5 
compounds. Parts (a) and (b) in Fig. 1-3 are three- 
dimensional plots of j4(k, to) versus w for fc-values in the 
ID Brillouin zone, whereas part (c) in all three figures 
summarize these results in the usual "band structure" ui 
versus k plot. The shaded areas represent regions with 
significant spectral weight, with the darker shading re- 
presenting more weight. Due to particle-hole symme- 
try, A(k = (k, )r),u) = A(k = (-7T — k, 0), —ui) where 
k = (k, ki_). In other words, one reflects k about k = 7r/2 
and ui about to get A(h,ui) for = 7r from A(h,ui) 
for ki_ = 0. This symmetry can be seen by comparing 
parts (a) and (b) of Fig. 1-3. Since this symmetry is not 
enforced by the Maximum Entropy procedure, it is only 
present approximately, and the extent to which the re- 
flected spectra at = do not match those at = 7r 
gives an indication of the accuracy of the analytic con- 
tinuation procedure. In the density plots in parts (c), 
we show only the = components. For t±/t = 2.0, 
there is a heavily weighted coherent band of width ~ 3t 
in the photoemission (to < 0) part of the spectrum. In 
the inverse photoemission spectrum (to > 0) there is ve- 
ry little total spectral weight. There are therefore two 
bands with significant spectral weight, one in the pho- 
toemission spectrum for = and one in the inverse 
photoemission spectrum for = 7r. Each band is about 
2t away from the Fermi surface, leading to a gap of ap- 
proximately At. Therefore, for large the system is a 
two-band insulator. 

In contrast, for t±/t = 1.0 and t±/t = 0.5, A(\t,ui) 
has substantial spectral weight in four bands, one in the 
photoemission spectrum (to < 0) for = and = 7r 
and one in the inverse photoemission spectrum (to > 0) 
for both values of In this regime, the system is thus 
a four-band insulator. For t±/t = 0.5 (Fig. 3) the spec- 
tral weight is present for both values of concentrated 
between k = and k = 7r/2 in the photoemission part 
and between k = tt/2 and k = 7r in the inverse pho- 
toemission part of the spectrum. For the isotropic case 
(ti_/t = 1.0, Fig. 2), there is some shift of spectral weight 
to the photoemission part for = and to the inverse 
photoemission part for = 7r. The maxima of the pho- 
toemission band in the = part occurs at k* = 7r for 
tj_/t = 2.0, at k* ra 0.7tt for t L /t = 1.0, and at k* = tt/2 
for ti_/t = 0.5. Therefore we would expect to see a ma- 
xima at k* 0.77T in a photoemission experiment on the 
isotropic ladder compounds, an experiment which has, to 
our knowledge, not yet been done. 



It is also important to note that the QMC A(k.,uj) 
results for t±/t = 0.5 as well as for the isotropic case, 
ti_/t = 1.0, display two general features which have re- 
cently been observed in both the ID 13 and 2D 14 cases 
using the improved Maximum Entropy techniques des- 
cribed above 16 . One feature is that A(h,u>) contains a 
rather dispersionless "incoherent background" extending 
over several t (~ 6t for U ft = 8 in 2D) in both the electro- 
nically occupied (to < 0) and unoccupied (to > 0) parts of 
the spectrum. The crucial structure, not resolved in pre- 
vious ID and 2D QMC simulations, is a dispersive struc- 
ture at low energies with a small width of the order of the 
exchange coupling J = At 2 /U . It is this latter coherent 
"band" which defines the gap A and which, for larger 
U (U/t ~ 10), is well-separated from the higher-energy 
background. As in the ID and 2D cases, the splitting in 
the low-energy band and the higher-energy background 
is especially pronounced near k = and k = 7r due to a 
relative weight shift from negative to positive energies as 
k moves through k* . 
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r 

FIG. 4. The spin-spin correlation function S(r) for t j_ = 2t, 
t± = t and t± = 0.5t calculated on a 2 x 32 lattice with 
U/t = 8 and open boundary conditions using the DMRG. 
The solid, dotted and dashed lines are asymptotic fits to the 
data using the simple exponential form Aexp( — r/{). 

B. Spin— Wave Theory 

At half-filling, the two chain Hubbard system at T = 
is a spin liquid with a spin gap and a finite spin-spin 
correlation length. The system is therefore less orde- 
red than either the ID system, which is at the critical 
point and has power law decay of the spin-spin correla- 
tion function and no spin gap, or the 2D system, which 
has long-range spin order and gapless spin-wave-like ex- 
citations. However, at large U jt the scale of the spin 
gap for the two chain system is set by Jj_ ~ At 2 j_/U, 



A 



and the spin-spin correlation length £ grows longer as 
the spin gap gets smaller. The range of the spin orde- 
ring can therefore be tuned by varying t±/t. This is 
illustrated by Fig. 4, in which we plot the spin-spin cor- 
relation function, S(r) = (— l) r (Mq A M r z A ) on a semilog 
scale for tj_/t = 2.0, tj_/t = 1.0 and t L /t = 0.5. Here 
M* x = («r,A,t — n r,\,i) is twice the z component of the 
on-site spin. These calculations were done using the Den- 
sity Matrix Renormalization Group method 8 (DMRG) 
on a system with open boundaries at zero temperature. 
For ti_/t = 0.5, the data was averaged over a number of 
different pairs of points for a given r in order to reduce os- 
cillations due to the open boundaries. For distances grea- 
ter than a few lattice spacings, the correlation functions 
are very well fit by a pure exponential form j4exp(— r/£), 
from which the spin-spin correlation length £ can be de- 
termined. For t L /t = 2.0, £ = 0.83, for t L /t = 1.0, 
£ = 4.3, and for t±/t = 0.5, £ = 9.5 lattice spacings. 
Therefore, the spins are almost completely uncorrelated 
along the chains for tj_/t = 2.0, but for tj_/t = 1.0 and 
for ti_/t = 0.5 are correlated on length scales of the order 
of the size of the 2 x 16 lattice studied here, although the 
long-range behavior is that of a gapped spin liquid in all 
three cases. 
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FIG. 5. The SDW gap calculated within the AFHF appro- 
ximation at T = 0, as a function of t±/t on a 2 x 100 grid for 
U/t = 2,4,8. 

It is therefore reasonable to calculate the single par- 
ticle spectral weight A(/k,u>) in two simple ways: first, 
when the spin gap is small and £ ps L, and we consider 
wavelengths smaller than £ and frequencies larger than 
the spin gap, it is useful to explore the consequences of 
a simple SDW approximation based on an ordered state, 
i.e. AFHF theory. Secondly, one can solve for the exact 
eigenstates on a rung and then consider states compo- 
sed of a product of exact rung singlets only and a delo- 
calized one-hole state in a background of rung singlets. 
(We call this the Local Rung Approximation, LRA.) The 
LRA starts with a state with uncorrelated rungs and then 
treats the interactions between rungs perturbatively, and 
is thus a good starting point when £ is small. As we shall 



see, the LRA gives a single particle spectral weight distri- 
bution that agrees well with QMC in the large t j_ regime, 
while AFHF qualitatively describes most features of the 
spectral weight distribution for small tj_. 

We first discuss the AFHF calculation in more detail. 
The Hartree-Fock one-particle Hamiltonian can be writ- 
ten as 

k,<7 

where ^ k J means that the sum runs only over the ma- 
gnetic zone (— 7r/2 < k < 7r/2, = 0, ir). The operators 
6£ a and a are given by the following transformation 18 : 

K,a = «kCk,<7 T fkCk+Q,<7, (6) 
K,a = ^kCk,<7 ± WkCk+Q,<7, (7) 

where Q = (tt, 7r), the upper (lower) sign corresponds to 
a =| (<t =|), and 

«k = ^(l + £ k /^k), (8) 

v k = ^(l-e k /E k ), (9) 
E k = VA 2 +£k 2 - (10) 

Here is the free particle energy (U = 0) defined in Eq. 
(2). The spin-density-wave gap, A, is self-consistently 
determined by the equation 

k 

The band structure of the AFHF Hamiltonian is given 
by i-Ek = ±\/ A 2 + £k 2 in the magnetic Brillouin zone. 
The SDW gap A, calculated by solving Eq. (11) self- 
consistently, is shown as a function of tj_/t for U/t = 2, 
4, and 8 in Fig. 5. The transition from a two-band metal 
to a band insulator in the U = system occurs at tj_/t = 
t_i_c = 2.0. As U is turned on for t < t± c in the AFHF, a 
gap develops in both the bonding and antibonding bands, 
leading to a four-band insulator. If t±_ is then increased, 
A will go to zero at approximately the noninteracting 
(U = 0) t ±c /t = 2.0, as seen for U/t = 2 in Fig. 5, 
leading to a transition from a four-band to a two-band 
insulator. For large U and small A ps U/2, as seen 
for U/t = 8 in Fig. 5, so the Coulomb splitting of the 
bands will dominate over U = band structure effects. 
In this case, the gap will go to zero, i.e. the system 
will undergo a transition to a band insulator only when 
the band structure splitting 2tj_, is of the order of the 
Coulomb splitting, so that t j_ c ps A ps U /2. 

The spectral weight within the AFHF approximation 

is 

A(k,u 1 ) = ul6(u 1 -E li ) + vl6(u 1 + E li ). (12) 
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We show the spectral weight A(k,u>) from AFHF for the 
same parameters as in Fig. 1 and Fig. 3 in Fig. 6(a) for 
t ± /t = 2.0 and Fig. 6(b) for t L /t = 0.5. The band split- 
ting 2A pa U in both cases, as can also be seen in Fig. 5. 
Since both the bonding and antibonding bands are split 
by this gap, there are two peaks in A(k,u>) as a function 
of w for each k, leading to a four-band insulator for both 
values of tj_/t. Since the band splitting is set by U rather 
than 2tj_, the spectral weight is almost evenly distribu- 
ted between two coherent bands for both tj_/t = 2.0 and 
ti_/t = 0.5. Therefore, for tj_/t = 2.0, as can be seen 
by comparing the QMC results in Fig. 1 with Fig. 6(a), 
the AFHF spectral weight distribution for tj_/t = 2.0 is 
completely wrong, with too much weight in the upper, 
<jj > 0, band for = and the u> < band for = ir. 
While the average positions of the AFHF bands, shown 
as thick solid lines in Fig. 1(c), are approximately the sa- 
me as those of the QMC bands, the band widths are too 
small by approximately a factor of two. For tj_/t = 2.0 
and U/t = 8, then, the AFHF calculation results in a 
four-band insulator in which there is long-range anti- 
ferromagnetic order, while the QMC results show that 
the system is a two-band insulator, with only very short 
range antiferromagnetic correlations. While the AFHF 
picture does produce a two-band insulator for t\_jt > 
4.5, the physical picture of the transition to this phase is 
quite different (see section C). 
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FIG. 6. Spectral weight calculated within the AFHF ap- 
proximation for (a) t±/t = 2.0 and (b) t±/t = 0.5 with 
U/t = 8, shown for k± = and k± = ir. 

For tj_/t = 0.5, as seen in Fig. 3 and Fig. 6(b), AFHF 
gives two dispersive bands in the = branch of width 
of order J ~ At 2 /U , similar to the coherent bands seen 
in the QMC data. However, in the AFHF, the single 
particle gap is somewhat too large, the weight distribu- 
tion extends too far towards k = 7r and too far towards 
k = for the photoemission and inverse photoemission 
parts of the spectrum, respectively, and the broad inco- 



herent background is not present. We have also carried 
out mean-field slave boson calculations 20 in order to try 
to improve on the AFHF picture. These calculations give 
a gap about 15% smaller than AFHF, and a bandwidth 
about 3% smaller, but do not qualitatively improve on 
the AFHF calculations. The same holds for tj_/t = 1.0, 
where the AFHF describes approximately the coherent 
part of the spectrum and also produces a too large single 
particle gap and a slightly different spectral weight distri- 
bution. Therefore, for tj_/t = 0.5 and for tj_/t = 1.0, AF- 
HF would give a reasonable description of the dispersion 
and general spectral weight distribution of the coherent 
spin-wave bands, if the gap were phenomenologically ad- 
justed to a smaller value, but fails to even qualitatively 
describe the spectral weight distribution for tj_/t = 2.0. 
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FIG. 7. Exact eigenstates for a two-site system represen- 
ting a disconnected rung, with U = 8 and t± = 2.0. The 
horizontal axis is the number of particles and the vertical 
axis shows the total energy ((H — fiN)), where the chemical 
potential fi = U/2 at half-filling. Each level is labeled with 
the total spin S and momentum k±. 



C. Local Rung Approximation 

For very large a better starting point is the limit 
of weak interaction between the rungs. In this limit, we 
split the Hamiltonian of Eq. (1) into H = Ho + Hi with 

Ho = -t± ^{c\ M c i 2a + h.c.) + n iM n iM , 



(13) 



where Ho denotes the non-interacting rung limit. The 
ground state of Ho is a product of individual rung ei- 
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genstates. Accordingly, we first diagonalize the Hubbard 
Hamiltonian on the two sites making up the rung in or- 
der to find the rung eigenstates. A schematic diagram of 
the exact eigenstates for a single rung is shown in Fig. 
7. The ground state for the half-filled rung i (two par- 
ticles per rung) is a spin singlet state with = and 

energy E /L = £™ ng = (u - ^V 2 + 16t\j /2 denoted 

by \Si). For large values of U, the energy of this two-site 
state is given by the exchange coupling — Jj_ ~ — At\/U . 
Note that \Si) contains terms which have double oc- 
cupied sites as well as the usual singlet construction 

(|U)-IU»/V2. 

We can then form an approximation to the ground 
state for the half-filled system by taking a state which is 
just a product of the rung states, 

\^o) = \S 1 )\S 2 )...\S L ). (14) 

For ti_/t = 2.0, the binding energy of a singlet formed 
on a rung is about four times lower than the energy of 
a singlet between neighboring sites along a chain, and 
therefore we expect |t/>o) to be a good approximation to 
the exact ground state. 

In order to obtain the u> < spectral weight for the 
half-filled system, we need to calculate the matrix ele- 
ment of ci^ct given in Eq. (3). In other words, we need 
matrix elements between a half-filled state and a state 
with one particle removed. We form the approximate one 
hole state by replacing one bond singlet state at rung I, 
\St), with the lowest energy one-particle state, the one 
with bonding (k±_ = 0) symmetry, \B^) (We remove a 
spin down electron for definiteness.) and define 

\£) = \S 1 )\S 2 )...\B n )...\S L ). (15) 

We delocalize the single-particle state with a plane wave 
ansatz by constructing the state 6,21 

\Mk)) = L-^j2 eM \ e )- ( 16 ) 

i=i 

Trivially, this state is an exact eigenstate of Ho- We will 
take this state as the ground state for the system with 
one hole with momentum k = (k,0). The spectral weight 
A(k,u>) at zero temperature (/? — ► oo) is then approxima- 
ted by using only the two states |t/>o) and |t/>i(fc)} in Eq. 
(3). For = 0, the energy dispersion u>(k) of A(k,w) 
for ijj < is given by 

uj(k) = {MHHo) - {Mk)\H\Mk)) - V 

= -\\j lP + m l +< -L -tAcosk, (17) 

w ith A = ( 1 + E 2 /2t ± f/(l + El/At\), E 2 = (U + 
\/U 2 + 16t jJ/2, and the corresponding spectral weight 
by |(i/>i(A;)|cfc ) jJi/>o}| 2 - The dispersion given by Eq. (17) 
for U/t = 8 and t L /t = 2.0 is plotted in Fig. 1(c) and 
the spectral weight for = is shown in Fig. 8. One 



can see that the position of the peak of the u> < LRA 
band, denoted LRA1 in Fig. 1(c), lays almost exactly on 
the QMC data. However, the position and the dispersi- 
on of the ijj > band, which were obtained in the same 
way as for u> < 0, are not exactly matched by the LRA1 
calculation. While this band is not very important in 
the sense that it contains very little spectral weight, we 
can understand how to improve the LRA1 calculation by 
considering the states of a four site cluster. 
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FIG. 8. Spectral weight using LRA1 calculation for the 
lo < and the LRA2 calculation for the lo > parts of the 
spectrum for k± = 0. Here t±/t = 2.0 and U/t = 8. 

In order to generate the inverse photoemission spec- 
trum, we need to calculate matrix elements between the 
half-filled state and a state with one additional particle. 
From Fig. 7, we see that a rung state with three partic- 
les and ki_ = has approximately twice the excitation 
energy of the = one-particle state which is relevant 
for the photoemission part of the spectrum. (Recall that 
due to the particle-hole symmetry at half-filling, this 
state will map to one in the inverse photoemission part 
of spectrum for k = 7r; one can see this symmetry in Fig. 
7.) Since the relevant three-particle single rung state is 
high in energy, configurations involving intrachain effects 
might also be important. We include the effect of such 
configurations by replacing two rung singlet states by the 
lowest energy state of five particles on four sites in the 
bonding channel. We can then form a delocalized plane 
wave state from this state as in Eq. (16). The results of 
this calculation are labeled as LRA2 in Fig. 1(c). The 
location and width of the = inverse photoemission 
band are more closely fit than by the LRA1 calculation. 
In Fig. 8 the spectral weight distribution of the LRA1 
calculation for the photoemission part of the spectrum 
(w < 0) and the distribution of the LRA2 calculation for 
the inverse photoemission part (u> > 0) is plotted. The 
result is in good accordance with the QMC data in Fig. 
1(a). 
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FIG. 9. Dynamic spin response function x s (q, w) for (a) 
= 2.0, gj. = 0, (b) = 2.0, gj. = ir, (c) tj./* = 1.0, 

g_i_ = 0, and (d) t±/t = 1.0, gj_ = ir. Again, [//i = 8 on a 
2x16 lattice. The weight of sizable structures is represented 
by the strength of shading in the shaded regions and the posi- 
tion of the maxima by the points with error bars. In (b), the 
solid line indicates the spin gap A_E™ ng of the two-site system 
and the dashed line the spin gap AEs for the 2 x 16 ladder 
at T = 0, calculated with a projector QMC algorithm. In (c) 
and (d), the solid lines are X^p A (q, q, w) calculated using the 
SDW approximation given by Eq. (19). 



III. SPIN AND CHARGE DYNAMIC 
CORRELATION FUNCTIONS 

In order to determine the nature of the low lying exci- 
tations, we also consider the spin and charge susceptibi- 
lities Xs,c(q,^) which, in a Lehmann representation, are 
defined as 



X S)C (q,c) = -^e-^(l-e-^)|(/|0 S)C (q)|/'}| 2 



(18) 



p+q, 



S+q.lS.l) aIld ^ 



with O s (q) = £p(< 

X2 P (j c p+q <j c p a- We calculate the two-particle dynamic 
response from the two-particle imaginary time Green's 
function, as in Eq. (4) and use the Maximum Entropy 
method described above for the analytical continuation 
of x SiC (q, w) to real frequencies. 

In Fig. 9 we show the dynamical spin-spin correlation 
function Xs,e(q,<^) for tj_/t = 2.0 and tj_/t = 1.0. For 
ti_/t = 2.0, the q±_ = component [Fig. 9(a)] has a broad, 
lightly weighted structure centered approximately at u> = 
3t, with the spectral weight vanishing for small q. For 
q±_ = 7r [Fig. 9(b)], there is a coherent, dispersive band 
with a width set by J k, At 2 /U = O.bt, a minimum at 
q = 7r and a maximum at q = 0. The minimum spin gap 
is approximately 0.8t, which agrees well with the value we 
obtain using Projector Quantum Monte Carlo, indicated 
by the dashed line on the plot. The spectral weight is 
most heavily concentrated around q = 7r. 

The dynamic spin response at half-filling measures the 
response of the system to spin triplet excitations. Since 
here t j_ is large, we consider the effect of triplet excita- 
tions on product rung states of the type considered in 
Eq. (14), in the LRA calculation. A triplet excitation 
will change the total spin of the state from S = to 
S = 1. Referring to Fig. 7, the only triplet excited state 
on a single rung has momentum 7r (corresponding to odd 
parity under chain interchange) leading to a change in 
momentum q±_ = 7r for the triplet excitation. The size 
of this triplet excitation AE 1 ™" 8 is marked on Fig. 9(b) 
by the solid line. This local triplet excited state can be 
moved to a near neighbor rung by a process which is se- 
cond order in Hi. The local triplet excited states can 
then be delocalized into a Bloch wave as in Eq. (15) and 
Eq. (16). For large U, one then obtains a dispersion re- 
lation of the form AE 1 ™" 8 + J cos q, where J ~ At 2 /U , a 
form similar to that obtained in Ref. 5 for the two chain 
Heisenberg model. The coherent band in Fig. 9(b) does 
have a minimum at q = 7r, and has approximately this 
form. 

Within the LRA picture, it is also possible to excite 
a q±_ = (even under chain interchange) triplet excita- 
tion via a "two-magnon" process, as discussed in Ref. 
5. This corresponds to making two local triplet excita- 
tions on rungs in the non-interacting rung picture, and 



8 



would lead to an excitation energy whose scale, in the 
Heisenberg limit, is set by u> ~ 2Jj_ rather than u> ~ Jj_. 
The overall position of the lightly weighted band seen for 
= in Fig. 9(a) is consistent with this two-magnon 
band, although the band is too lightly weighted and our 
resolution too low to extract a dispersion relation. 

For the physical relevant isotropic case (tj_/t = 1.0), 
shown in Fig. 9(c) and (d), the spin response looks quite 
different. (Because the tj_/t = 0.5 results are qualita- 
tively similar to the tj_/t =1.0 results, we show only the 
isotropic (tj_/t = 1.0) case here.) In Fig. 9(c) and (d), 
there are dispersive bands in both the q±_ = and q±_ = 7r 
branches, with the position of the peak going to a finite 
minimum at q = (0,7r) and q = (tt, 0), and appears to 
vanish as q — ► (tt, 7r). However, as shown in Fig. 4, the 
correlation length is 4.3 lattice spacings in this regime, 
and there should be a spin gap. Using DMRG calculati- 
ons on lattices of 2 x 8 to 2 x 32 sites 8 , we estimate the 
spin gap to be of order 0.12t. Due to finite size effects 
and finite resolution in the analytic continuation, this 
gap is too small to resolve in Fig. 9(d). Near q = (0, 0) 
[Fig. 9(c)], the dispersion of the peak is hard to discern 
because there is very little spectral weight. This lack of 
spectral weight at q = (0,0) is present in all the dyna- 
mic charge and spin correlations and is due to a selection 
rule that comes about because S)C (q = 0) in Eq. (18) 
commutes with the Hamiltonian, leading to a vanishing 
of the matrix element of S)C (q) as q — ► 0. Therefore, 
the system has relatively long-range spin order in this 
regime and, to within the resolution of our calculations, 
shows the characteristics of an ordered state with gapless 
excitations. Thus it is interesting to compare the QMC 
spectra to the results of spin-wave theory calculations. 

In order to extract the low-lying spin excitations using 
spin-wave theory it is necessary to consider the RPA 
transverse spin susceptibility Xrpa^' k', ui). It is obtai- 
ned by applying the random-phase approximation to the 
response function Xo~ (k, k', ui), which is calculated di- 
rectly in real time using the AFHF ground state 
calculated in section II B: 

Xj-(q,q',*) = ^<n|TS+(*)S: q ,(0)|fi>, (19) 

with = S x ± iSy . Due to the broken spin rotational 
symmetry of Xrpa( < 1' q> w ) con t ams a gapless mode, 
as predicted by the Goldstone theorem. Following the 
procedure of Schrieffer et al. 18 , one obtains: 



XR PA (q>q'> w ) = 



(20) 



where [1 — Uxo (P;q'; w )] 1 is a matrix inverse of a 2 x 2 
matrix in q-space and 

Xo""(q, q',w) = <5(q- q')xo~~(q> w ) 

+%-q' + Q)x$"(q,w), (21) 



with Xo~ (q> w ) an d Xq (q; w ) given by the usual "bub- 
ble" diagrams and printed in detail in Ref. 16 [up to a 
misprint of a factor of 2 in XQ - (q ; w )]- The dispersion 
of XRPA(qi q; w ) i s indicated as a solid line in Fig. 9(c) 
and 9(d). The RPA spin-wave dispersion goes to zero at 
q = (0, 0) and q = (tt, 7r), and goes to a finite minimum 
at q = (0,7r) and q = (tt, 0), consistent with the QMC 
data. The spin-wave velocity, the width of the bands, 
and the gaps at q = (0,7r) and q = (tt, 0) are also in 
reasonable agreement with the QMC data. 



Ul/t 




FIG. 10. The dynamic charge response function x c (q, w) 
for the same parameters as in Fig. 7. In (a), t±/t = 2.0 
and q± = ir, the solid line represents the charge excitation 
energy on a single rung A_E™ ng , and the line with solid dots 
represents the minimum excitation energy constructed from 
the single particle spectral weight in Fig. 1. In (b) and (c), 
t±/t = 1.0, and q± = and 7r, respectively. The lines with 
solid dots represent the minimum excitation energy of the 
RPA charge susceptibility Xrpa(1j 1j ^ j calculated with the 
method of Ref. 16. 
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The QMC results for the dynamic charge correlation 
function are shown in Fig. 10 for tj_/t = 2.0 and tj_/t = 
1.0. For ti_/t = 2.0, almost all of the spectral weight 
is in the q±_ = 7r component, so we do not show the 
q±_ = component. There exist two important features: 
one is a heavily weighted, relatively flat band at u> ~ lOt 
with heaviest weight near q = 0. This band becomes 
somewhat incoherent as q increases. The second is a flat, 
dispersive, less heavily weighted band with a minimum 
of order 5t near q = ir. The spectral weight in this band 
extends from about q = tt/A to q = ir, and the size of the 
charge gap is set by this lightly weighted lower band. 

In the large t±_ limit, one can understand the structure 
of the charge response from the LRA picture described 
in section II C. In the single rung picture, a charge exci- 
tation will occur through a transition to an excited state 
conserving the number of particles on the rung and the 
total spin; in other words, the important transition will 
be in the middle column of Fig. 7, from the low-lying 
5 = state to the higher S = states. There are two 
possible charge excited states on the rung, one with mo- 
mentum = (even parity) and one with momentum 
= 7r (odd parity). In a single rung picture, an optical 
transition from the = ground state to a = 
excited state is forbidden because the q±_ = density 
operator O c (0) commutes with the Hamiltonian. This 
selection rule forbidding a q±_ = optical transition re- 
mains present when the rung charge excited states are 
constructed as in Eq. (15) and delocalized in a state like 
that in Eq. (16). This is why there is almost no spectral 
weight in the q±_ = portion of the charge response for 
ti_/t = 2.0. Of course, the system is not exactly in a LRA 
state, so there will be some higher order processes that 
will introduce a very small amount of spectral weight into 
the q±_ = branch. The energy of the q±_ = 7r single rung 
transition, indicated by a solid line on the Fig. 10(a), gi- 
ves an excitation energy that agrees well with energy of 
the heavily weighted region at q = 0. 

In order to qualitatively understand the origin of the 
dispersion of the charge response, one can consider the 
possible particle-hole excitations within the one-particle 
band structure given by A(k,u>) in Fig. 1. There will be 
significant amplitude in the two-particle charge response 
when there is significant amplitude for a transition at a 
particular q = (q,q_i) and u> for a particle-hole excita- 
tion built up from the one-particle spectral weight. For 
example, to understand the heavily weighted amplitude 
at q = (0,7r), one has to integrate over all transitions 
from the photoemission band in Fig. 1(a) to the inverse 
photoemission band in Fig. 1(b) which transfer this mo- 
mentum. Since the single-particle bands are sharp and 
parallel, one should obtain a single well-defined peak for 
q = (0, 7r), which we see in Fig. 10(a) by considering exci- 
tations between the u> < 0, = band, and the u> > 0, 
= 7r band. In addition, the transition at q = (0,7r) 
has odd parity between the chains, and is thus allowed by 
the selection rules for the density operator as q — ► 0. As 
the parallel component q is increased, one can see that 



there will be a continuum of excitation energies which 
gets wider as q increases. The minimum excitation ener- 
gy (position of the lower band) from this construction is 
shown in Fig. 10(a) by the line with solid dots. Howe- 
ver, the excitation energy obtained is consistently smaller 
than the energy of the lowest heavily weighted band from 
QMC. Using the single particle bands to construct the 
two-particle excitations is equivalent to calculating the 
charge response using the lowest order "bubble" diagram, 
but with exact single-particle propagators, neglecting all 
particle-hole interactions. The particle-hole interactions 
on the rung, which are included in the rung eigenenergies, 
thus raise the charge excitation energy by a substantial 
amount 22 . 

For ti_/t = 1.0, the charge response looks quite diffe- 
rent. As shown in Fig. 10(b) and (c), there is substan- 
tial spectral weight for both q±_ = and q±_ = 7r. For 
q —> (0, 0), the density operator has even parity, causing 
the matrix element and thus the spectral weight to va- 
nish, whereas at q — ► (0, ir), the density operator has odd 
parity so that optical transitions are allowed and there is 
spectral weight. In both channels, the size of the charge 
gap is approximately At, and there is a broad structure 
of width ~ It. For q±_ = 0, most of the spectral weight 
occurs as a dispersive peak whose energy increases with 
increasing q, whereas for q±_ = 7r two peaks seem to con- 
tribute to the spectral weight distribution. The peaks 
are not well-defined enough over a range of q to extract a 
dispersion, but the upper peak is heavily weighted near 
q = 7r, at ijj ?s 9t. 

We can qualitatively understand the broad incoherent 
structure of the charge response by considering particle- 
hole excitation in the single particle A(k,u>) in Fig. 2. 
There are four dispersive bands and a broad background, 
so there should be weight in both the q±_ = and q±_ = 
7r branches of the charge response, and broad structure 
above a minimum excitation energy, which we see in Fig. 
10(c) and (d). From the single particle bands in Fig. 3, 
one can estimate the minimum particle-hole excitation 
energy to be ~ At. In Fig. 10(c) and (d), the spectral 
weight near this minimum excitation energy is suppressed 
due to the particle-hole vertex. We have also carried out 
a calculation of charge response Xrpa^i l' > w ) within the 
antiferromagnetic RPA approximation described above, 
using the SDW dispersion in Eq. (10), and also find a 
relatively broad structure above the charge gap for both 
gj_ = and q±_ = ir. We plot the minimum excitation 

energy of Xrpa^' l' > w ) m Fig. 10(b) anc ^ ( c ) as lmes 
with solid dots. This line is located in the middle of the 
broad band in both plots. 

The spin and charge response functions for tj_/t = 0.5 
show the same general features as in the isotropic (t_i/t = 
1.0) case and are therefore not shown here. The spin 
susceptibility x s (q, w) is also qualitatively identical to 
the RPA result Xrpa( < 1' *l> w ) wl th a smaller spin velocity 
than in the tj_/t =1.0 case. The charge susceptibility 
Xc(q,w) shows a clear dispersive band centered around 
the low-lying RPA excitations in the q±_ = channel, 
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whose energy increases with increasing q. In the q±_ = 7r 
channel, there is a broad structure of width ~ 8t with 
again two peaks in the spectral weight distribution. 
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IV. CONCLUSION 

In summary, the single and two-particle dynamical 
properties of the two chain Hubbard model at half-filling 
can be understood by starting from two limits: the limit 
of non-interacting rungs treated exactly, which gives a 
good starting point for the large t±_ case for which the 
spin-spin correlation length along the chains is less than 
a lattice spacing, and the limit of an antiferromagneti- 
cally ordered state, which gives a good starting point for 
the small t j_ case when the spin-spin correlation length is 
large. The dynamical properties in the two regimes look 
quite different. For the large t±_ regime, the remnants of 
the level transitions of the two site system representing a 
single rung, suitably broadened into bands, can explain 
the major features of the single-particle, spin and char- 
ge responses. In the small t±_ case, calculations based 
on an antiferromagnetically ordered starting point, such 
as antiferromagnetic Hartree-Fock theory and spin-wave 
theory give a good qualitative picture of the coherent 
spin-wave part of the single particle spectral weight and 
of the two-particle spin dynamic correlation function. In 
addition, there is a broad incoherent band in the single 
particle spectral weight similar to that found in recent 
numerical work on the ID and 2D systems. 

We also have shown the single particle spectral weight 
and the spin and charge response functions for the physi- 
cal relevant, isotropic case (tj_/t = 1.0). The results are 
qualitatively similar to these in the small t j_ region and 
can therefore be understood from an antiferromagnetic 
Hartree-Fock or spin-wave theory picture. The real lad- 
der compounds have approximately the same coupling 
strength parallel to and perpendicular to the chains, and 
it will therefore be interesting to compare these isotropic 
ti_/t =1.0 results with future experiments. 
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